Super-shell structures and pairing in ultracold trapped Fermi gases 
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We calculate level densities and pairing gaps for an ultracold dilute gas of fermionic atoms in 
harmonic traps under the influence of mean field and anharmonic quartic trap potentials. Super- 
shell structures, which were found in Hartree-Fock calculations, are calculated analytically within 
periodic orbit theory as well as from WKB calculations. For attractive interactions, the underlying 
level densities are crucial for pairing and super-shell structures in gaps are predicted. 
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Ultracold atomic gases have recently been used to cre- 
ate novel quantum many-body systems such as strongly 
interacting high temperature superfluids of fermions, 
Bose-Einstein condensates, Mott insulators in optical lat- 
tices, etc. These lab phenomena have a strong over- 
lap with condensed matter [lj, nuclear [2] and neutron 
star physics [3]. Finite fermion systems such as atoms in 
traps, nuclei, helium and metal clusters, semiconductor 
quantum dots, superconducting grains, etc., have addi- 
tional interesting quantum structures such as level spec- 
tra, densities and pairing. These will be observable as 
temperatures are further lowered in atomic trap experi- 
ments. The high degree of control over physical param- 
eters, including interaction strength and density, makes 
the atomic traps marvelous model systems for general 
quantum phenomena. 

The purpose here is to calculate the level spectra, 
densities and pairing for zero-temperature Fermi gases 
in harmonic oscillator (HO) traps with anharmonic and 
mean field perturbations, and to show that novel super- 
shell structures appear in both level densities and pairing. 
In calculating level spectra by analytical periodic orbit 
theory and WKB as well as numerical Hartree-Fock, we 
also relate these different theoretical approaches to one 
another. 

We treat a gas of N fermionic atoms of mass m in a HO 
potential at zero temperature, interacting via a two-body 
interaction with s-wave scattering length a. We shall 
mainly discuss a spherically symmetric trap and a dilute 
gas (i.e. where the density p obeys the condition p|a| 3 <C 
1) of particles with two spin states of equal population. 
The Hamiltonian is then given by 



H = 



N 

E 



2m 



mw 2 r- + U(ri) 
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We will consider both external anharmonic potentials of 
the form U — er 4 and particle interactions: Ufa) — 
(2-Kh 2 a/m) 3 3 ( r i~ r j)- When interactions are weak, 
the latter can be approximated by the mean field poten- 
tial 

,_, , 2irti 2 a , . 

U(r) = p(r) . (2) 



For a large number of particles and U — the Fermi 
energy is Ep — npfuo where tif — Uf — 3/2 ~ (3N) 1 / 3 
is the HO quantum number at the Fermi surface. The 
HO shells are highly degenerate with states having an- 
gular momenta I = np, rip — 2, mod(riF, 2), due to the 
U(3) symmetry of the 3D spherically symmetric HO po- 
tential. However, interactions split this degeneracy. In 
the Thomas-Fermi (TF) approximation (see, e.g., [4|]) the 
Fermi energy is 
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The density p{r) = k F (r)/3ir 2 is 

p(r) = PQ (l-r 2 /R 2 TF f 2 



(3) 



(4) 



inside the cloud r < Rtf = a osc y/2nF, where pa = 
(2nf) 3 / 2 /3n 2 a 3 sc is the central density [5||. For conve- 
nience we set the oscillator length a osc = ^Jh/muj = 1 in 
the following. 

Taylor expanding the density and thereby also the 
mean field of Eq. J2]) around the center gives 



p{r) ~ p 1 



TF 



l /R 



TF 



(5) 



the first term will simply incorporate a constant shift 
in energies whereas the term quadratic in radius renor- 
malizes the HO frequency as io e g = ujy/l — 6napo/R FF . 
The third term is quartic in radius and is therefore also 
of the same form as the external potential 



U(r) 



er 



(6) 



with e — (3irh 2 a/ Am) po / R TF . Both the pure quartic 
potential and the mean field potential of Eq. ([2]) are an- 
harmonic and change the level density by splitting the I 
degeneracy of the HO shell nj? at the Fermi surface. 

We will now calculate analytically the level spectra 
from perturbative periodic orbit theory for the quartic 
potential and subsequently within semiclassical WKB 
wavefunctions for both the quartic and the mean field 
potential of Eq. We will start with repulsive inter- 
actions where pairing is not present. 
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In periodic orbit theory [6J], the level density can be 
written (to leading order in fi. -1 ) in terms of a perturba- 
tive HO trace formula 0, S| 



9(E) 



E 2 



(Hlj)- 



1 + Re (-l) k Me l27rk 



kE/hu 



(7) 



For the unperturbed HO (U = 0) the modulation factor is 
M. = 1. For a quartic perturbed potential, as in Eq. ([6]), 
the modulation factor was calculated in [8j 



M = 



f e -i2ka /h-i-rr jl , e -i3ka /h+in /2\ 

ka V / ' 



(8) 



with a — ettE 2 /h 2 uj 3 , being a small classical action. The 
two terms arise from the change in actions for the circle 
and diameter orbits respectively due to the quartic po- 
tential Q. The resulting level density can be written in 
the factorised form [§] 



9(E) = 



E 2 



(hwf 
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Here, the first term is the average level density, the 
cosine factor gives the rapid HO shell oscillations (mod- 
ified by the perturbation) which, however, are slowly 
modulated by the sine factor resulting in a beating pat- 
tern. Moreover, the non-perturbed HO limit, equivalent 
to M. = 1 in Eq. (7]), is recovered in the limit of |e| — > 0, 
where the U (3) symmetry is restored. The k = 1 term 
in Eq. J9]) gives the major oscillations in the level density 
and is shown in Fig. 1 (a) . The beating pattern or super- 
shells is clearly observed. The shell oscillations vanish 
when the argument of the sine in Eq. |(9|) is an integer 
S = 1,2,3,... times tt, i.e. \e\E 2 /2(hw) 3 = S. This gives 
the supernode condition 



n F = E/fkj = y/2Shu/\e\ 



(10) 



We now turn to an alternative calculation of the level 
density with WKB. The splitting of the HO shells degen- 
erate levels I — np, np — 2, mod (n F , 2) in the shell np 
by the mean-field potential can be calculated perturba- 
tively in the dilute limit. An excellent approximation for 
the radial HO wave function with angular momentum I 
and (n F — l)/2 radial nodes in the HO shell when np ^> 1 



is the WKB one 



Hn F i(r) 



2 sm(ki{r)r + 6) 



.1/2 



(r)r 



(11) 



between turning points r\ = hp ±^h 2 F -l{l + l). Here 
hp = rip + 3/2 and the WKB wave number k[(r) is 



kf(r) =2h F -r 2 -1(1 + 1) /r 2 



(12) 



When np ^> 1 the wave function has many nodes 1 <C 
I <C np and the oscillations in 72.^ (r) can be averaged 
(sin 2 (fc;(r)r)) = 1/2 The phase 6 is then unim- 

portant. The single-particle energies for the anharmonic 
potential of Eq. ([6]) are simply 



E„ 



hpfnu 
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U(r)\K nFl (r)\ 2 r 2 dr (13) 



= / + -^-\dr = £ - [3h 2 F - 1(1 + 1)] .(14) 
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It is special for the quartic perturbation that the level 
energies are linear m 1(1 + 1). The resulting level spacing 
increases as (2Z + 1) just as the level degeneracy for 5*0(3) 
symmetry. Therefore the level density is constant within 
the bandwidth 



D — |-En F ,2=0 ^ E n 



en 2 F /2 



(15) 



on energy scales larger than 2D /np but smaller than D. 
The level density vanishes between the bandwidths of 
two neighbouring n shells and therefore it generally has 
a strong oscillatory behavior as shown in Fig. Q] (a). Its 
amplitude is largest when D ~ huj/2. However, when 
D ~ hw the level density is constant and the oscilla- 
tory behavior vanishes. This phenomenon repeats when 
D = Shw since the level spectra then overlap 5* times. 
With the bandwidth of Eq. lfl~5j) under this condition, we 
obtain exactly the same supernode condition as for pe- 
riodic orbit theory, Eq. (jTUJ) . We conclude that Craig's 
perturbative periodic orbit theory 0| is in exact agree- 
ment with perturbative WKB for a quartically perturbed 
spherical symmetric HO in three dimensions. 

We now turn to the slightly more complicated mean 
field potential of Eq. |(2|). Its level spectrum can also be 
calculated from the WKB wave functions of Eq. (fTTjl . 
Inserting them in Eq. lfT3|) . we obtain 



E nf 



hphijj = 2/ (37r 2 



ah 3 J 2 huj I . 



Here, the integral I is 
/ = 



(\ - x^l-l{l + \)/h 



3/2 



dx 



(16) 



(17) 



where x = (r 2 —hp) / \Jh? F — 1(1 + 1). This integral is / = 
tt for / ~ np and / = 8\/2/3 for I = 0. The bandwidth is 
therefore 



D = 2/ (3tt 2 ) anj 2 ha; (&V2/3 - tt 



(18) 



Inserting this bandwidth in the supernode condition D = 
Shw gives 



(19) 



(8%/2/3 - tt) 2/ (3tt 2 ) anf = S . 



For example in the case 2na = 1 the supernodes S = 
1,2,3,.. should occur when np ~ 28,44,58, etc. The 
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Hartree-Fock (HF) calculations of the oscillating part of 
the total energy, which is proportional to the level den- 
sity at the Fermi level @, result in slightly higher su- 
pernodes, as in Fig 1 (b). The differences arise because 
the WKB calculations are perturbative in the interaction 
strength, whereas in the HF calculation the MF poten- 
tial U includes a large scattering length which, e.g., leads 
to corrections for the effective oscillator frequency. Also 
for the purely quartic term the perturbative approach 
underestimates the exact supernodes (see Fig. 3 of 0). 
For weaker interactions 2ita = 0.1, the first supernode 
S = 1 should occur at np = 130 according to the condi- 
tion of Eq. (fl9|) , in closer agreement with the HF result 
of Fig. Ql(c). 

For comparison, the Taylor expansion of the mean field 
potential leads to the supernode condition of Eq. lfT0|) 
with e = (3Trh 2 a/4m)p /RT F - It differs from Eq. JUJ 
by the prefactor, which is ~ 34% smaller. It is a better 
approximation to expand e.g. around r = Rtf/2^/2 = 
v / tTf/2, where the corresponding prefactor is only ~ 8% 
smaller, such that the supernode in Fig 1 (c) is predicted 
to np — 137. Now expanding I of Eq. (fl7|) for small 
I <C tif, one finds 



(8/3)V2- l 2 /V2n 



F • 



resulting in the level spectrum [1( 



E 



,,,,7-fificd = ■^ran 3 J 2 hu! 



16 1(1 + 1) 



(20) 



•(21) 



This level density is constant at low I as for the potential 
in Eq. (fT4|) . However, near I ~ np the density of lev- 
els is slightly smaller as can be seen from the bandwidth 
corresponding to Eq. {2TJ, which is ~ 12% larger, for a 
given np, than the bandwidth of Eq. lfl8|) . That the level 
density is not completely constant within the bandwidth 
has the effect that a small periodicity remains even at 
the super-shell condition D — Shw. Therefore the shell 
oscillations do not disappear completely at the supern- 
odes, as can be seen in Fig.[T](b,c), whereas for the purely 
quartic case (a) the oscillations disappear completely at 
the supernodes. 

Most atomic traps are not spherical but cigar shaped 
(prolate) with u z ^,ui±. The unperturbed HO energies 
E = n z huj z + n±hai± will generally lead to a constant 
level density for energy scales larger than huj z but smaller 
than Huj± . When the oscillator frequency ratio cu± /lu z is a 
rational number, level degeneracies and larger oscillations 
will occur on the scale Huj z . Interactions will, however, 
smear this level density. In any case, super-shell structure 
is not expected as in the spherical symmetric case. In 
very oblate traps w z ^> u>± the mean field potential is 
effectively two-dimensional and quadratic, i.e. it does 
not split the HO shells [13, E|- Thus we may expect 
strong oscillations in the level density on the scale hai±, 
but again no super-shell structure. 
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Figure 1: (color online) The upper figure (a) shows the leading 
term (k = 1) of the oscillating part of the perturbative level 
density of Eq. (J9} as a function of uf = E/hu>, for the case of 
an external potential V = V H o+er 4 with e = 2/40 2 « 0.0013. 
The middle and lower figures (b,c) show the oscillating part 
of the total energy according to a numerical HF calculation 
with interaction strength 2%a — 1 and 2%a — 0.1, as a 
function of the HO shell number (h = u — 1). This illustrates 
qualitatively that a supernode, e.g. at «f = 40, can be due 
to interaction (b) and/or an additional quartic term to the 
HO potential (a). 



Attractive interactions lead to pairing by an amount 
that is exponentially sensitive to the underlying level den- 
sity near the Fermi surface (3, EH EH, 0] ■ The level den- 
sity is the same for repulsive and attractive interactions 
except that the levels are reversed when the sign of e 
(Eqs. dSJ) and fl3|) or a is changed (Eq. ([16|). Therefore 
we can use the level densities and bandwidths calculated 
above for pairing calculations. Pairing in finite systems 
is described by the Bogoliubov-de Gennes (BdG) equa- 
tions [lH and take place between time-reversed states. 
As shown in [3] these states can be approximated by HO 
wave functions in dilute HO traps as long as the gap does 
not exceed the oscillator energy, A^hui. Solving BdG for 
such finite systems is numerically complicated and we 
shall therefore apply further simplifying approximations, 
namely that the pairing gap A n i and the wavefunction 
overlap matrix elements vary slowly with level I in a shell 
n. Both approximations are fair for the trapped atoms 
as argued in [ll|| and deviations can be understood. As 
result we arrive at a much simplified gap equation 



1 = 



G_ 

^2" 



-2n f 



9(E) dE 



y/(E-»)* + A(»y 



(22) 



Here, the supergap G = Z2\f2nF\a\fko /15ir 2 was cal- 
culated in [ifjj as the pairing gap when all states in 
a shell can pair; this is the case for a region of in- 
teraction strengths and particle number where the gap 
is large as compared to the level splitting, yet small 
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stronger at midshell than at endshell, 
fewer states to pair 



30 40 
n =(3N) 1 ' 3 



Figure 2: (color online) Multi-shell pairing gaps for a HO trap 
with an additional quartic term in the potential with e = 
2/40 2 , i.e. for the level density of Fig. 1 (a) with supernodes 
at riF — 40, 40\/2 « 57, etc. The interaction strength is a = 
—0.05 (top red curve), a = —0.03 (middle blue curve, with the 
inset figure around the first supernode) and a = —0.01 (lower 
green curve). In the inset plot it is clearly seen that the local 
minima for I ~ tif and I ~ before the supernode turns 
into local maxima after the supernode, as a consequence of 
overlapping shells. The dashed (red) line is the multi-shell gap 
A = G/(l — 2Gln(nF)/hu>) for a = —0.05 and the upper/lower 
thin solid line (black) are the single mid-/end-shell pairing for 
a — —0.01 (see text). 



compared to the shell splitting ttirf. A(/x) = A„; is 
the gap at the Fermi surface. g(E) = n F /D is 
the level density within each bandgap D around ev- 
ery shell n — 0, 1, ...,~2njr but vanishes between 
the bandgaps. The gap equation thus reduces to 
1 = (G/D) I D dE/^{E + nhu- fx) 2 + A 2 . The 

chemical potential fi can be determined from the level 
spectrum; as we gradually fill particles into the shell np 
at the Fermi surface, /i increases from nphui to npftio+D . 
The cut-off n^2n F in the sum of the gap equation models 
as a first approximation the more rigorous regularization 
procedure described in Ref. that is required for a 
delta-function pseudo-potential. 

By solving this simplified gap equation of Eq. (|22| . 
we find that it still contains and displays the essential 
interplay between the variation in level density and pair- 
ing. To illustrate the super-shell structure in pairing, 
we take the strongly anharmonic trap potential used for 
the level spectra in Fig. 1 (a), and calculate the pairing 
arising from a weak attractive scattering length a < 0. 
For sufficiently weak interactions such that pairing only 
takes place in the shell at the Fermi surface, we obtain 
the expected result from the gap equation: A = G when 
D < A, whereas for D > A we get A = D exp(-D/2G) 
midshell (// = n F hw + D/2) and A = 2Dexp(-D/G) 
endshell (/i = npfnu or [i = npfnu + D). Pairing is thus 



where there are 
and strong shell oscillations fol- 
low as shown in Fig. 2. For stronger interactions, pairing 
also takes place between states in shells around the Fermi 
shell and Eq. (22]) gives: A = G/ (1 - 2G\n{n F ) /hu) for 
small bandwidth [14|. In Fig. 2 this curve is compared 
with the finite bandwidth result, which has strong oscil- 
lations except at the supernodes where the level density 
is continuous. At a supernode D = fuv and the gap equa- 
leads to a gap A = 2n_pfiwexp(— hui/2G) 



tion (122|) 

In summary, level densities, shell-oscillations and 
super-shell structures in anharmonic traps calculated 
from numerical Hartree-Fock and analytical periodic or- 
bit theory as well as WKB were found to match to lead- 
ing order. Analogous super-shell structures were found 
in pairing from an approximated BdG calculation. The 
mean field in atomic nuclei also have a large anharmonic 
potential and the HO shells start to overlap (the first 
supernode) already for heavy nuclei with np ~ 5 — 6. 
The interplay of level spectra and multishell pairing is, 
however, difficult to disentangle in nuclear pairing due to 
strong spin-orbit effect and small particle number. Ul- 
tracold atomic traps, however, provide ideal systems for 
observing the rich quantum structures such as level den- 
sities and pairing. 

Discussions with Matthias Brack on periodic orbit the- 
ory, Ben Mottelson on (nuclear) shell theory and pairing, 
and proof reading by Joel Corney, are gratefully acknowl- 
edged. 
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